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Nowadays Nuclear Magnetic Resonance diffusion (dNMR) measurements of water molecules in 
heterogeneous systems have broad applications in material science, biophysics and medicine. Up to now, 
microstructural rearrangement in media has been experimentally investigated by studying the diffusion 
coefficient (D(t)) behavior in the tortuosity limit. However, this method is not able to describe structural 
disorder and transitions in complex systems. Here we show that, according to the continuous time random 
walk framework, the dNMR measurable parameter a, quantifying the anomalous regime of D(t), provides a 
quantitative characterization of structural disorder and structural transition in heterogeneous systems. To 
demonstrate this, we compare a measurements obtained in random packed monodisperse micro-spheres 
with Molecular Dynamics simulations of disordered porous media and 3D Monte Carlo simulation of 
particles diffusion in these kind of systems. Experimental results agree well with simulations that correlate 
the most used parameters and functions characterizing the disorder in porous media. 



I n the last years, much attention has been devoted to develop techniques for investigating biological tissues in 
I vivo by using non invasive procedures. Among them, Nuclear Magnetic Resonance diffusion (dNMR) has 
I proved to be a successful tool with its ability to provide detailed information about material 1 " 7 and tissue 8 ' 9 
properties. The measurement of molecular diffusion is an effective experimental way of probing biological and 
porous material structures, since the geometrical complexity of heterogeneous systems can be characterized by 
the way in which small molecules (typically water) diffuse within them. 

A complex porous structure implies a multi-scale hindering of diffusing molecules so that its study and 
characterization can reveal the fingerprint of the medium geometry. 

In conventional dNMR studies, structural complexity in heterogeneous media has been investigated by means 
of tortuosity indices, e.g. the obstruction factor 0 D W . However, we show here that the parameter 0 D is not suitable 
for detecting the most relevant information on global structural complexity (disorder) and structural transitions. 
In this work we focus on an innovative aspect of dNMR that allows, in the hindered diffusion regime, to describe 
the disorder properties of a sphere packed medium, to monitor its structural modifications and, if present, to 
reveal structural transitions. 

In particular, we show that the exponent a, a measurable dNMR parameter quantifying the anomalous regime 
of the time- dependent diffusion coefficient D(t), is related to well known statistical quantities characterizing the 
disorder in porous media (the bond-orientational order parameters 11 , and the pair correlation function). Notably, 
a is able to detect structural transitions and to classify different kinds of disorder. 

To achieve this goal, we compare the result of a dNMR experiment made on a sample of random packing 
monodisperse micro-spheres with the outcome of numerical 3D simulations resulting from an original "collage" 
of validated algorithms to yield essential results in the characterization of complex systems 12 " 14 . 

Specifically, to highlight the physical meaning and relevance of a we show here a detailed description of 
structural disorder in random packing of mono-dispersed spheres performing both generative Molecular 
Dynamics simulations of disordered porous media and extensive 3D Monte Carlo studies of particles diffusion 
in these kind of systems, eventually connecting the numerical results to experimental dNMR measurements. 

We explain the results of particle diffusion in the framework of the continuous time random walk (CTRW), 
which has been demonstrated to be a general and effective theoretical tool able to quantify anomalous diffusion in 
both laboratory- and field-scale systems 1516 . 

The diffusion process of solute particles in dilute homogeneous and isotropic solutions is described by the 
celebrated Einstein -Smoluchowski equation: 
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(Ar(t) 2 )-(Ar(t)) 2 = 2d Dt 



(I) 



where the expression on the left of Eq.(l) is the mean square dis- 
placement (MSD) of solute molecules describing the spatial fluctua- 
tions around their coherent motion, d is the effective topological 
dimension of the medium in which the transport process occurs, 
and D is the solute diffusion coefficient. 

In highly heterogeneous media the migration of solute particles is 
affected by various combined mechanisms effectively determining 
trapping into geometrically restricted zones. This leads to a broad 
spectrum of effective "local" solute particle stopping times, for which 
the relation (1) is no more valid. In general, natural systems show a 
range of heterogeneity length scales 151718 , £, between a minimum £ min 
and a maximum £ max . Both theoretical and experimental studies 18-23 
have shown that three distinct diffusion regimes of solute particles 
can be identified in heterogeneous media. For £<£ m [ n , diffusion is 
not sensible to heterogeneities and the relation (1) holds with a 
constant D equal to the bulk diffusivity D 0 of the medium. 

On the other hand, when t>£ maXi caging and trapping effects due 
to heterogeneities are massively averaged out and relation (1) holds 
again but with an effective constant diffusion coefficient <D 0 , 
which takes into account the average effect of traps and barriers 24 . 

Finally, on intermediate scales £ min < £ < £ max the diffusion is 
sensible to heterogeneity variations and may show a more complex 
behavior, known as non-Fickian or anomalous diffusion, which in the 
case of statistical isotropy takes the form of MSD =2d D(t)t, with a 
D(t) that can be often well fitted by a power law 151718 



(Ar(t) 2 )-(Ar(t)) 2 = 2dTf 



(2) 



The parameter a is defined as the anomalous diffusion exponent and 
r is a generalized diffusion coefficient. Usually the term anomalous 
diffusion refers to an asymptotically (in time and scale) anomalous 
regime, i.e. the t a behavior should persist for t—>™. This behavior can 
result from an infinite hierarchy of heterogeneity length scales. 
However, many complex real porous materials, such as those inves- 
tigated in this study, exhibit a finite hierarchy of heterogeneity length 
scales which is a sufficient condition for anomalous diffusion at short 
times crossing over to normal diffusion at long times. Therefore, in 
this work, instead of defining away the transient anomalous 
diffusion, we analyze and quantify the initial period of anomalous 
diffusion. 

This scenario introduces a generalized effective D(t) that in the 
case of Eq. (2) takes the simple form: 



D(t) = Tf- 



(3) 



When 0 < a < 1 the diffusion process is sub -diffusive, for a > 1 is 
super- diffusive and when a = 1 the diffusion is normal and the 
relation (1) is recovered. 

Real heterogeneous porous media are expected to host ordinary or 
sub-diffusive processes (i.e. a < 1). The transition from the inter- 
mediate anomalous diffusion regime to the asymptotic ordinary 
regime is characterized by the crossover time f *, related to the largest 
heterogeneity length scale l max ~ y/t*D^. 

Conventional dNMR structural characterization of heterogeneous 
media is based on the measure of the so called obstruction factor 0 D 
= DJD 0 , i.e. the inverse of the tortuosity index of the medium 1012 . In 
general, 0 D depends primarily on the shape and orientation of the 
geometrical restrictions as well as on the packing ratio (j) 25 . The exact 
analytical derivation of the 0 D dependence on <j> is a tricky math- 
ematical problem. However, empirical relations with adjustable 
parameters have been introduced 25 " 27 in analogy to well known 
theoretical results about simple cases. 

We classify the degree of disorder in heterogeneous monodisperse 
sphere packing systems at different angular scales through the widely 
used 11 ' 28 " 30 bond-orientational order parameters Q/. Their values are 



invariant under rotations of the reference coordinate system and take 
characteristic values at the crystalline ordered phase. Therefore they 
can be used to quantify the type and degree of rotational symmetry in 
the system. 

The local bond-orientational order parameter is defined by 11 



,(*) = 



1/2 



(4) 



with Yi >m spherical harmonics, 0 jk and <pj k polar angles of the vector 
between particles j and k, and (. . .)^ enn m denoting the average over 
the nearest neighbors j of k. 

The global bond-orientational order parameters, Q b are often used 
to quantify global structural disorder 30 . In order to compute Q/ the 
average (. . .)^ ewn m in equation (4) is substituted by the average over 
all pairs of nearest neighbors. 

We restrict our analysis to / = 4 and 1 = 6, which bear particular 
importance in the case of cubic and icosahedral symmetries, respect- 
ively 11 . Since the first non vanishing Qj for the icosahedral, hexagonal 
and cubic symmetries is Q 6 , we chose it as an indicator of the degree 
of local order in the system 28 " 30 . We normalize Q 6 to its value in the 

perfect fee crystal and defined the quantity Q = as the normal- 
ized bond-orientational order parameter. In the infinite volume limit, 
Q ranges from 0 (complete disorder) to 1 (perfect fee order) and 
therefore can be used to quantify the disorder degree. Hence, a suit- 
able definition of nearest neighbor in a disordered medium is needed. 
The medium we considered is composed by monodisperse hard 
spheres. We define a nearest neighbor of a given sphere k any other 
sphere with its center closer than a certain distance r min to the center 
of k. We chose r min as the length scale of the first relative minimum of 
the pair correlation function g 2 (r) defined as usual as 



#(r) = 



2nNr 2 p 0 \2"'*; 



1 



r,) 



(5) 



with AT being number of spheres, p 0 = N/V, Vthe total volume of the 
system and (...) the ensemble average. 

The expression p 0 g 2 (r)dV delivers the probability to find the cen- 
ter of a particle in the volume element dVata distance r around the 
center of a given sphere. 

Details on the numerical estimation of g 2 (r) are given in the 
Supplementary materials (SM). 

Results 

dNMR and numerical experiments. We investigated a system made 
up of hard and impenetrable monodisperse spherical micro -beads in 
water at different sphere packing, as a proxy for a simple porous 
medium. 

dNMR experiments were performed on five samples made of 
beads of nominal mean diameters: 30.0, 20.0, 15.0, 10.00 and 
6.00 urn mono-dispersed in water solution. The measured values 
of the packing ranges from <j> = 0.35 to <j> = 0.56 (see Methods for 
details). 

The systems used in the NMR experiments were numerically 
reproduced by 3D Molecular Dynamics simulations of N = 500 
mono-sized quasi-hard spheres with periodic boundary conditions 
on the horizontal plane and reflective conditions on the bottom of the 
vertical direction. The numerical simulations yielded systems with a 
broad range of packing values 0.35 < <j> < 0.74. The face centered 
cubic (fee) crystal, corresponding to a perfectly ordered close packed 
situation with </> = 0.740, was built with the corresponding primitive 
vectors. A schematic representation of some of the systems simulated 
in this study is shown in Fig. 1. 

Finally, the diffusion of water molecules in the interstitials of such 
systems, was modeled by a Monte Carlo 3D random-walk process 
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with reflecting boundary conditions at the surface of the spheres that 
constitute the inaccessible medium. 

Geometry-based analysis of numerical results. The average 
structure seen by a generic particle of the system described by g 2 (r) 
displayed in Fig. 2, shows a full agreement with the predicted 
theoretical regimes found in literature 28 ' 31 . In all cases, we observe a 
pronounced peak at a distance r = a, with a the sphere diameter that 
corresponds to the distance of the nearest neighbors in contact. For r 
> a, the probability to find neighbors decreases reaching a minimum 
around r ~ 1.4a. In fluid-like systems, theoretically for </>< 
0.55 28 ' 31,32 , the g 2 (r) is known to oscillate with decreasing ampli- 
tude. By increasing the sphere packing and staying far from a 
crystallization regime, one gets a random jamming scenario, theo- 
retically for 0.55 < (j) < 0.64 28 ' 31 ' 32 , and the g 2 (r) displays two peaks at 
r~\f3a and r ~ 2a, i.e. at the distances of the second and third 
neighbors of a closed-packed crystalline structure, still oscillating 
for larger values of r. The value (j) ~ 0.64 defines the maximum 
random jamming (MRJ) configuration 31,32 . Finally, in systems with 
larger sphere packing (j) > 0.64, signs of crystallization become more 
evident and the g 2 (r) displays multiple narrow peaks that are the 
signature of a crystal structure on distances larger than the first 
neighbor distance 33 . The graphs in Fig. 2 clearly indicate that in 
our numerical samples a jamming transition from an unjammed 
liquid-phase to a disordered jammed phase takes place around the 
MRJ packing. 

A deeper analysis of the local packing configurations can be 
obtained by exploring spherical local neighborhoods and can be 
useful to identify cooperative effects and possible structural transi- 
tions. From the study of the geometric nearest neighbors number 
distribution (GNNd) and the local spheres density fa distribution 
(LSDd) we found (see SM) a decreasing number of realizable con- 
figurations when the system approaches the random close packing. 
This behavior can be related to the onset of cooperative effects and of 
a possible phase transition from unjammed liquid to disordered 
jammed phase. These results are in good agreement with recent 
measurements of granular systems 34 " 37 , supporting the existence of 
a transition from random fluid system to glass or amorphous solid, 
geometrically frustrated on local crystal structure, but disordered on 
large scales (polycrystal) 38 . 

The average number n c of spheres in contact with any given 
sphere, is able to estimate more precisely the two different sphere- 
density thresholds associated with the two phase transitions known 
in literature 39 " 41 , occurring in our simulated disordered systems. The 
first transition, from fluid to an intermediate phase that is rigid but 
unstressed 39 , is identified by the rigidity percolation threshold (RPT) 



we observe at <j) ~ 0.52 (the theoretical value 39 " 41 is found at <j> ~ 0.54). 
The second transition, from this marginal state to a geometrically 
frustrated solid state, is identified by the stress percolation threshold 
(SPT) that we observe at <j> ~ 0.62 (the theoretical value 39 " 41 is the 
MRJ packing found at (j) ~ 0.64) (see SM). 

Geometrical frustration, polycrystals and effective anomalous 
diffusion. It has been often argued that the driving mechanism 
generating amorphous structures could be the geometrical frustra- 
tion 42 . In order to test whether such a frustration mechanism occurs 
in our simulated systems in the amorphous solid phase for (j) > 0.62, 
we searched for configurations that were locally close packed at 
densities larger than the crystalline ones (details in SM). The geome- 
trical frustration would have an important role in amorphous 
systems if a non-negligible fraction of such configurations would 
be found. We found that the geometrical frustration plays a negli- 
gible role when (j) < 0.70, where the fraction of geometrically 
frustrated particles f GF was found to be less than 10%, in agree- 
ment with previous experimental works 40 ' 41 . Conversely, at high 
values of sphere packing (j) > 0.70, in the amorphous disordered 
phase, local sphere arrangements with high local densities, such as 
the icosahedral packing, play an important role. In that case we found 
40% <f GF < 60% (see SM Fig. S5). This means that in our highly 
packed and disordered hard sphere systems an abrupt transition 
from an unjammed liquid-phase to a disordered jammed phase 
occurs and frustration arises while the continuous translational sym- 
metry breaks locally. In this case, due to the creation of frustrated 
areas between locally ordered configurations, systems at high (j) 
values exhibit a broad spectrum of typical heterogeneity length 
scales (see Fig. 3). This is in fact confirmed by the size distribution 
of the interstitial voids that has a variance a\ d showing a minimum at 
(/) ~ 0.62 and increasing monotonically (see SM). 

The appearance of multiple heterogeneity length scales by frustra- 
tion suggests, according to the CTRW framework, that in very dense 
amorphous packing an effective anomalous diffusion process should 
be observed. In particular, on the basis of the arguments given above, 
we expect to measure an a < 1 at these (j) values. The behavior of D(t), 
D vs diffusion time f, for all the simulated systems, i.e. fee ordered 
(black points) and disordered (blue points), is reported in Fig. 4, 
where we notice that the diffusion in disordered systems exhibits a 
dynamics with a more anomalous character than in ordered ones. 
Since the decrease of D(t) in ordered systems is due to hindering 
effects by construction, the corresponding more accentuated 
decreasing behavior in the disordered regime may be ascribed to 
the complexity of the system itself. In particular, the decreasing rate 
at high (j) values is more pronounced in disordered samples, while the 
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Figure 1 | Hard sphere systems. The spherical regions used for the calculation of some of the statistical quantities described in the text is shown for 
different hard sphere packing ratios and degrees of disorder. The intersection between the spherical surface and the spheres is depicted in blue color. 
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Figure 2 | Pair correlation function. The g 2 (r) is shown for selected values of hard sphere packing. The horizontal axis 
sphere diameter. Three different scenarios can be distinguished: i) a fluid-like behavior for (j) < 0.55; ii) random jamming 
amorphous solid for (/) > 0.64. 
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width of the transient anomalous diffusion region, marked by the red 
straight lines in Fig. 4, is broader. Moreover, Doo gets lower values at 
higher sphere densities, with a trend more noticeable in disordered 
systems than in ordered ones. All these findings confirm the theor- 
etical and numerical results of previous works 12 dealing with the 
obstruction factor 0 D . 

Local bond orientational order parameters analysis. In order to 
establish the nature of the local internal organization emerging at 
high packing values, and to measure how close the packing is to 
an ideal crystal structure, we determined the bond orientational 
order parameters Q 4 and Q 6 for all considered systems. Below 
(j) < 0.62 we found (see SM) less than 5% of local configurations 
with symmetry compatible with either bcc (body centered cubic) 
or fee and less than 10% with local symmetry compatible with hep 
(hexagonal close-packed). As (j) increases, the local close-packed 
configurations hep and fee become more and more significant, 
reaching more than 30% of the total configurations as (j) 
approaches 0.70. This, together with the previous discussion on 
the transition of/ GF , suggests that a local crystallization process is 
setting up around (j) ~ 0.70. To highlight the onset of this liquid- 
polycrystal transition we plot in Fig. 5- (a) the behavior of Q vs (j), 
which shows that Q is very sensitive to any kind of crystallization 
transition and increases significantly when order appears around 
(j) ~ 0.70. In particular, Q increases monotonically with (j) as the 
system approaches the rigidity percolation threshold (RPT) at (/) ~ 
0.52, indicating that the increase in density yields an increase of 
short-range order in the fluid. A discontinuous jump in Q is 
observed at (j) ~ 0.62. At this value of 0 a significant amount of 



structural order appears, in perfect agreement with the stress 
percolation threshold (SPT). Hence, we deduce that the simu- 
lated systems undergo an abrupt transition at (j) ~ 0.62, from 
an unjammed liquid-phase to a disordered jammed phase. 




Figure 3 | Artistic representation of Polycrystal structure. Globally 
disordered structures comprised of frustrated grains locally ordered in 
crystalline configurations. They are characterized by a wide distribution of 
voids characteristic length scales that reflect disorder properties. This 
suggests, according to the CTRW framework, that an effective anomalous 
diffusion process should be detected in such systems. This means that a 
straight relation between the properties of structural disorder and 
diffusion in the hindered regime exists. Unlike 0 Dy the dNMR parameter a 
embodies information about the properties of structural disorder in an 
effective way and allows to identify the structural transition from rigid but 
unstressed random jamming to stressed polycrystal structure (see Fig. 5-b 
and Fig. 6). 



SCIENTIFIC REPORTS | 3 : 2631 | DOI: 1 0.1 038/srep0263 1 



4 



0.3 




log[t] (simulation units) 



Figure 4 | Time-dependent diffusion coefficient. Results of Monte Carlo 
simulations for the set of fee ordered (black points) and disordered (blue 
points) hard sphere packing at increasing packing ratio (j) ranging from 
0.33 to 0.74 (top to bottom). Data corresponding to the ordered systems 
are suitably shifted on the vertical axis for readability (D 0 is the same for all 
these systems). Red lines are power-law fits and quantify the scaling 
behavior of D(t) in the anomalous diffusion transient. Fee ordered hard 
sphere arrangements at packing ratio less than its maximum value reached 
at a closed packed geometry, i.e. for 0 < 0.74, were obtained by decreasing 
the diameter a of the spheres and letting their centers lie on the fee lattice 
points. The three diffusion regimes are clearly distinguishable for both 
ordered and disordered systems: free diffusion for t < 10 in simulation 
units (a time simulation unit correspond to ~2/is); transient anomalous 
diffusion when 10 3 ;$ t S 10 4 in the case of ordered systems, and 
10 2 ^ t ^ 10 4 in case of disordered systems; hindered diffusion in 
tortuosity limit for t> 10 5 . 

Discussion 

The main result of our study is depicted in Fig. 5 -(b) (a vs <j> graph). 
The a parameter is able to identify the structural transition occurring 
at <j) ~ 0.62 (fluid - amorphous solid transition identified by the 
abrupt change of slope) in similar way to the order parameter Q. 
Therefore, the diffusion of solute particles is sensitive to the disorder 
structural transition and can in principle be used as an effective 
experimental probe to monitor it. By looking the a experimental 



results displayed as yellow diamonds compared with simulations, 
it is possible to distinguish the random fluid phase (a > 0.96), the 
mixed fluid-solid phase (0.94 < a < 0.96) and eventually the 
amorphous solid state (a < 0.94). 

Conversely, 0 D vs <j> showed in Fig. 6, because of its high sensitivity 
to local properties of the medium, is not able to characterize any 
transition, and its use is limited to the discrimination between fully 
ordered and disordered systems 25 " 27 . 

As a first application of the method exposed here, we clarify the 
nature of the micro-beads systems experimentally studied by dNMR 
techniques. The measured values of indicate that these samples 
were in a fluid-like and fluid- solid like phase since 0.35 < <j> < 0.55. 
The measured values of 0 D ensure that the samples were in a dis- 
ordered state, well reproduced by our simulations. Moreover, the 
measured values of a show that the samples lay far from the jamming 
transition point. Finally, we can assert that three of the five samples 
were characterized by a random fluid-like microscopic state and the 
remaining two were in a mixed phase of mixed fluid and solid. 

In conclusion, we highlighted here important properties of the a 
parameter that can be easily measured by using dNMR experiments 
with faster procedure than in conventional D(t) behavior studies 10 . 
Thanks to the ability of NMR to be non-invasive and non-destruct- 
ive, a measurements could lead to the identification and monitoring 
of changes in tissue pathology or material structural-properties. In 
order to demonstrate that anomalous diffusion of particles in het- 
erogeneous porous media can provide information about their struc- 
tural complexity, a procedure composed of generative Molecular 
Dynamics simulations and 3D Monte Carlo studies was described. 
These simulations highlight the relation between the most suitable 
parameters and functions used to describe disorder in heterogeneous 
porous systems. Numerical and experimental results show an excel- 
lent agreement, confirming the reliability of a to characterize 
experimentally the structural disorder, opening the way to future 
developments and applications involving structural transition in het- 
erogeneous media. 

Methods 

Samples preparation and experimental set-up. Five 10 mm NMR tubes were filled 
up to a volume of approximately 2 cm 3 with polystyrene beads (Microbeads AS, 
Norway) of nominal average diameters of 30.0, 20.0, 15.0, 10.00 and 6.00 um, mono- 
dispersed in a solution of polyoxyethyle-sorbitan-mono-laurat (Tween 20) at 10" 6 M 
and deionized water (conductivity ~ 10" 6 Q"Vcm), and investigated during 4 months 
after their preparation at fixed temperature of 293 K. The measured values of the 
packing ranged from (j) = 0.35 after two weeks, to 0 = 0.56 at the end of the 
observation after four months. 




Figure 5 | (a-b) Structural order phase diagrams, (a) Normalized bond-orientational order Q vs packing ratio 0, for all the simulated disordered packing 
(green squares) as well as for the fee crystal (red cyrcle). (b) anomalous diffusion exponent a vs packing ratio (/). Values of a inferred by Monte Carlo 
simulations of solute diffusion in the fee ordered systems (red circles) and in the disordered systems (green squares) reported in Fig. 4, are showed 
together with dNMR experimental values (yellow diamonds). dNMR experiments were carried on five samples with 0.33 < 0 < 0.56. Experimental error 
bars for the dNMR data were estimated by a fitting procedure employing a Levenberg-Marquardt algorithm. Values of a were determined from the red 
line fits showed in Fig. 4 by using the definition of Eq. (3). In both (a) and (b), the planes are subdivided into colored regions, associated to the structural 
phases of the hard sphere systems. The rigidity percolation threshold (RPT) and the stress percolation threshold (SPT) are highlighted. 
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Figure 6 | Experimental dNMR O d values as a function of (f) compared 
with simulations results. Values of the obstruction factor O d inferred by 
Monte Carlo simulations of solute diffusion in the fee ordered systems (red 
circles) and in the disordered systems (green squares). Values were 
determined from the asymptotic behavior of D( t) at long diffusion times. 
Yellow diamonds represent values experimentally determined by dNMR, 
as described in the Methods. The three colored curves represent the fit of 
the semi-empirical relations derived by Lerman 26 , Boudreau 27 and 
Rayleigh 25 with the values got from simulations. In particular, the Lerman's 
relation is O d = [A(l - 0)^ m ]~", with A = 1.00 ± 0.13, m = 1.28 ± 0.11 
and n= 1.00 ± 0.10 for the ordered systems and A = 1.00 ± 0.15, m = 1.32 
±0.14 and n = 1.19 ± 0.13 for the disordered ones; the Boudreau's 

relation is Od = - — — A — j^, with C 1 



l-Clog(l-0 
ordered systems and C = 0.467 

Rayleigh's relation is Od 



0.3266 ± 0.0053 for the 



0.011 for the disordered ones; and the 
1 30 
1-0 V 2 + a^ + b^° 
0.10, b= 0.32 ± 0.11 and c= 0.20 ± 0.10 for the ordered systems and a = 



, with a = 0.53 : 



0.64 ± 0.12, b = 0.17 ± 0.11 and c = 0.54 ±0.11 for the disordered ones. 



The effective fractional exponent a was measured by collecting the NMR pulse field 
gradient (PFG) signal attenuation, S{q, A), as a function of the diffusion time A and by 
using the asymptotic expression of the Fourier transform of the motion propagator 
for the sub-diffusive regime derived in ref. 43: 



S(q,A)cce 



(6) 



which holds when g 2 <C ^ is kept constant at varying A, with K a generalized 

diffusion constant and q the wavevector defined by q = — yg3, where y is the gyro- 

2n 

magnetic ratio of the water proton, g is the pulse magnetic field gradient strength and 
3 is the pulse magnetic field gradient duration. A spectroscopic PFG Stimulated Echo 
(PFG-STE) sequence with 5 = 2.2 ms, g = 0.10 T/m (i.e. ^—11240 m" 1 ) along x, y 
and z directions, repetition time TR = 2.5 s, number of average NS = 32 and 48 values 
of A in the range (10-500)ms was used to extract the Gi i=x ^ z values. Then the mean 
value, a, was computed by averaging over the three directions. A second spectroscopic 
PFG- Stimulated echo experiment was used to obtain experimental 0 D values. A/8 = 
400/2.2 ms, TR = 2.5 s, NS = 16 and 48 gradient amplitude steps from 0.026 to 
1.02 T/m along x, y and z directions were used to measure the effective diffusion 

coefficient at long time, Doo, from the slope of In a t low q values along 

|_S(q = 0,A)J 

each direction and then averaging over the three directions. A Spin Echo imaging 
sequence with TR = 2.5 s, 48 values of echo time, TE, in the range (3-130)ms, NS = 
16, slice thickness STH = 1 mm, field of view FOV = 8 X 8 mm and an in plane 
resolution of 62.5 urn was used to estimate 0 in each sample, as described in ref. 44. 

Molecular Dynamic simulation of 3D hard-sphere packing. We simulated 
numerically and investigated two different sorts of obstructed media: i) systems made 
of 500 mono-sized hard spheres ordered on a fee lattice, at different values of sphere 
packing, ^e{0.330, 0.340, 0.380, 0.420, 0.470, 0.520, 0.570, 0.620, 0.680, 0.740}; 
ii) systems made of 500 mono -sized hard spheres randomly displaced, at different 
values of sphere packing, 0e{O.35O, 0.382, 0.460, 0.465, 0.538, 0.572, 0.573, 0.575, 
0.614, 0.627, 0.635, 0.640, 0.670, 0.701, 0.724, 0.737}. In the case of the fee lattices, 
the different packing were obtained by changing the sphere diameters while leaving 



the sphere centers on the fee lattice points. The disordered hard sphere systems were 
obtained by a Molecular Dynamics (MD) approach. The hamiltonian of the system 
was taken as 



(7) 



where m is the mass of particle i, Ui the external potential energy and U 2 the pair 
potential. The external field was modeled by a constant gravitational field whose 
corresponding force acted along the vertical axis z. In order to simulate interactions 
between hard spheres, the pair potential was chosen harmonic with a suitably high 
harmonic constant. The corresponding force acts along the radial direction when the 
distance between two spheres centers is lower than the sum of their radii: 



-ViI7 2 = 



■ K(r-d)-®{r- 



-d) 



(8) 



with d = q { — qj the distance between the centers of spheres i and j, r the sum of their 
radii (their diameter in this case of identical spheres), K ~ 10 4 simulation units the 
harmonic constant and @ is the Heaviside step function. In order to avoid unrealistic 
spheres bumpings, a dissipative force, F y , was introduced to have a system of 
harmonic oscillators in overdamped regime: 



F y = y{vd) j2 0{r-d) 



(9) 



where v is the difference between spheres velocities and y is the damping constant, 
chosen equal to 0.9 in simulation units. 

Moreover, in order for the system to reach a static final state of equilibrium, a 
viscous friction force, proportional to the velocities of spheres, was also introduced 
and applied after a proper thermalization time (to obtain a disordered system the 
viscous friction force has to be applied quite early). 

The MD simulations were performed with periodic boundary conditions on the 
horizontal x, y plane and reflective conditions on the bottom of the simulation box. 
The initial configuration was chosen from a random poissonian configuration of the 
spheres with null velocities. We obtained our final systems of spheres with different 
disorder degree by: 1) changing the density of the initial random configuration, 2) 
pinning a subset of spheres near the bottom, 3) changing the "switch on" time of 
viscous friction force. Simulation procedure terminated when the total kinetic energy 
of the system went below a cutoff value of 10 ~ 8 in simulation units. 

3D random flight simulation. Diffusion in 3D obstructed media, coinciding with the 
network of pores of the hard sphere systems, was modeled by random walkers moving 
among obstructions as a random flight, where each step vector 3r is uniformly 
distributed on the surface of a sphere of fixed radius |<5r|. 

Once the displacement has been extracted, if the final position of the considered 
particle is empty the displacement is accepted, otherwise the particle remains in its 
initial position, and another particle is randomly chosen to move. In other words, 
attempted steps that would lead the particle into obstructed regions, i.e. inside one of 
the packed hard spheres, are rejected. This procedure simulates the reflecting 
boundary condition at the obstructions surface 12 . Every simulation runs from 10 4 to 
10 6 time steps and every run was repeated 100 times with a different initial disposition 
of 500 independent and non-interacting walkers for run. Thus, the ensemble average 
(...) was computed by averaging over 5 X 10 4 random-flight trajectories with random 
initial positions. 

In order to avoid any artifacts in simulation results, the length of the step vector 3r 
is chosen |<5r| ^ a/100, where a is the smallest dimension of the obstructing objects, 
i.e. the diameter of the hard spheres. This is due to the fact that, although at low values 
of 0 diffusion characteristic parameters as 0 D are virtually independent of 1 3r\ , at high 
value of 0 (e.g. (j) > 0.70) 0 D was found to increase linearly 25 with decreasing 1 3r\ . The 
length of \3r\ chosen is sufficient to avoid this kind of simulation artifacts. 
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